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. In this paper, we investigate the solutions of a generalized fractional diffusion equation that 

extends some known diffusion equations by taking a spatial time-dependent diffusion coefficient 
and ./V-dimensional case into account, which subjects to natural boundaries and the general initial 



condition. In our analysis, the presence of external force is also taken into account. We obtain 
explicit analytical expressions for the probability distribution and study the relation between our 
solutions and those obtained within the maximum entropy principle by using the Tsallis entropy. 
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> 1 I. INTRODUCTION 

■ Recently, the diffusion equations that generalize the usual one have received considerable attention due 

to the broadness of their physical applications, in particular, to the anomalous diffusion. In fact, fractional 
diffusion equations and the nonlinear fractional diffusion equations have been successfully applied to several 
physical situations such as percolation of gases through porous media [1], thin saturated regions in porous 
media [2] , standard solid-on-solid model for surface growth [3] , thin liquid films spreading under gravity [4] , 
in the transport of fluid in porous media and in viscous fingering [5] , modeling of non-Markvian dynamical 
processes in protein folding [6], relaxation to equilibrium in system (such as polymer chains and membranes) 
with long temporal memory [7] , and anomalous transport in disordered systems [8] , diffusion on fractals [9] , 
and the multi-physical transport in porous media, such as electroosmosis[f 0-f 1]. The properties concerning 
these equations have also been investigated. For instance, in [12] boundary values problems for fractional 
diffusion equations are studied, in [13] a fractional Fokker-Planck equation is derived from a generalized 
master equation, in [14] the behavior of fractional diffusion at the origin is analyzed and a connection between 
the Fox H functions and the fractional diffusion equations was investigated in [15]. Also a generalization 
of Brownian motion to multidimensional anomalous diffusion is considered by using fractional differential 
equation in [16], and a fractional radial diffusion was considered in [17]. 

In this direction, we dedicate this work to investigate a fractional diffusion equation which employs 
space and time fractional derivatives by taking a time-dependent diffusion coefficient, an external force and 
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./V-dimensional into account. More precisely, we focus our attention on the following equation 

_ pM = ^ ^^_ { ^-i2, (x>t _^)__|y, (a . )tr} _ ___ {3 M-i F{x)pMh (i) 

where < 7 < l,p,v G R, Af is the number of dimensions, and the diffusion coefficient is given by 
D(x,t) = D(t)\x\~ e , which is a spatial time-dependent diffusion coefficient, and F(x) = — dV Q^ is an 
external force (drift) associated with the potential V(x). Here, we use the Caputo operator for the fractional 
derivative with respect to time t and the Riesz-Weyl operator for the fractional derivative with respect to 
spatial x [18]. We work with the positive spatial variable x. Later on, we will extend the results to the entire 
real x-axis by the use of symmetry (in other words, we are working with d/d\x\ and d^ -1 / d\x\^~ ). Also, we 
employ, in general, the initial condition p(x, 0) = p{x) (p(x) is a given function), and the boundary condition 
p(x — >±oo, t) — > 0. For Eq.(l), one can prove that f^™dxx Ar ~ 1 p(x,t) is time independent (hence, if p(x,t) 
is normalized at t = 0,it will remain so forever). Indeed, if we write Eq.(l) as d] = — l/x J ^~ 1 d x (x J ^~ 1 J), 
and, for simplicity, assume the boundary conditions J(±oo) = 0, it can be shown that f^^dxx-^~ 1 p(x,t) 
is a constant of motion (see [19] and references therein). Note that Eq.(l) recovers the usual radial diffusion 
equation with memory effect for 7 = 1. When (p, 7, v) = (2, 1, 1) and Af = 1, Eq.(l) recovers the standard 
Fokkcr-Planck equation in the presence of a drift taking memory effects and A/"-dimensional into account. 
The particular case F(x) = (no drift) and D(x,t) = DS(t) with (/i,7) = (2,1) has been considered by 
spohn [3]. The case D(x,t) = Dd{t)\x\~ e with (p,v) = (2,1) and the case D(x,t) = D(t) with (p,v) = (2,1) 
have been investigated in [20] and [21], respectively. In [22], the solution for Eq.(l) with the boundary 
condition p(0,t) = p(L,t) = was investigated. 

Explicit solutions play an important role in analyzing physical situations, since they contain, in principle, 
precise information about the system. In particular, they can be used as a useful guide to control the 
accuracy of numerical solutions. For these reasons, we dedicated to this work to investigate the solutions 
to Eq.(l). We consider different scenarios involving the diffusion coefficient and external force. Firstly, in 
the case of the absence of external force, we consider a spatial time-dependent diffusion coefficient given by 
D(x,t) = l?|a;| _6 't a ~ 1 /r(a) and D(x,t) = D\x\~ B 8{t). Secondly, the presence of external force is investigated, 
which is given by F(x) oc xlxl™" 1 . Thirdly, a mixing between the spatial and time fractional derivatives 
case is investigated. Lastly, we consider a particular case as 7 = 1, D(t) = D5(t), F(x) = —ICx, and 6, p, v 
arbitraries. In all the above situations, Eq.(l) satisfies the initial condition p(x,0) = p(x) (p{x) is a given 
function), and the boundary condition p(±oo,t) = O.Thc remainder of this paper goes as follow. In Sec. 2, 
we obtain the exact solutions for the previous cases. In Sec. 3, we present our conclusions. 



II. EXACT SOLUTIONS FOR DIFFERENT CASE 



In this section, we start our discussion by considering Eq.(l) in the absence of external force with 
D(t) = Dt^/Tia) ( D{t) = DS(t)),(p,u) = (2, 1) and 7, 9 arbitrary. For this case, Eq.(l) reads 

^p(x,t)=^dt'D(t-t')J^^x^-^ x p(x, t ) } . (2) 
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Here, we use the Caputo operator [18] for the fractional derivative with respect to time t. By employing the 
Laplace transform in Eq.(2), we obtain 



d 







__D( s) __ {a ^-i-^ p {X: s)} _ s i~ p ( x , S ) = -s-y-^ix, 0), 



(3) 



where p(x, s) = J£{p(x, t)}, D(s) = J/? {£>(£)}, and Jzf {/(£)} = f Q dte st f(t) denotes the Laplace transform 
of the function /. This equation can be solved by Green function method [23]. By substituting 



p(x,s) = I dx'x' N - 1 g{x-x',s)p{x') 



into Eq.(3), which yields 



d 



—D^^-ix^-^^-Gix, a)} - a ig(x, s) = -s^Gix, 0) 



d = 



dx 



dx 



In order to solve Eq.(5), it is convenient to perform the transform [24] 

y = A{s)x\ g(x lS )=y 5 Z{y) 
to translate Eq.(5) into the second-order Bessel equation as 



,d 2 Z dZ 



If 



IS 



S+l/v 



Oy* + y aH - (A + V )Z{V) = - sAW-W ^AM™ 



with parameter A 2 under the following conditions: 

v = = -[- 

2 W v l D(s) 1 



1 , 2 + 8-N . 2 + 8-N 

2 A = d = 

2 + 9 ' 2+9 ' 



(4) 



(5) 



(0) 



(7) 



(8) 



Since Eq.(5) should fit the boundary condition Q(±oo,t) = 0, i.e. Q(±oo,s) = 0, we get the solution of 
Eq.(5) 



g(x,s) = C(s)y s K x (y), 



(9) 



where K n {x) is the modified Bessel function of second kind, and C(s) can be determined by the normalization 
of G(x,t), i.e. /°° dxx^-igfa s) = ±. After some calculations, we obtain 



nl \ 2 + 8 1 S 1 i 2+e+M 2 + e-M 2 , S 7 1 2+£ 
g(x,Sj = u ( (- ) 2 ) 2 + 6 \x\ 2 K i+e-M ( (- )z\x\ 2 ). 

2+a v 2 + 9 D(s) 



r(^2 + e^ D (s 

here, we have used the formula 



' dy • y«K,(ay) = 2^a~^ r(l±^)r(i±^; 



(10) 



(11) 



Case 1. D(t) = D6(t), i.e. D(s) = D. 



Since K\(x) = \H^\ 



— |(-A/2,1)(A/2,1) 



, we can get the Laplace inverse of g(x,s) by applying the 



property of the Laplace inverse of Fox function [25], which yields 

1 



0M = J^( 



2T{^y{2 + 8YD^ 1 1 ' 2 



{2 + 8) 2 DV 



(1-0.7) 

(0,1), (^#,1) 



(12) 
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FIG. 1: The behavior of p(x,t) in Eq.(13) is illustrated by considering C(^,0,M)Q{x,t) versus ^ 2 +e)' 1 Dt~i ^ or typical 
values of 7 and 9, where C(~{,6,N) = 2r ^/ + ( 2+fl » ((2 + 9) 2 DV) M /{ - 2+e) . Here, for simplicity, we consider p(x,0) = 



where ff™/ 1 



(ai,Ai),...,(a p , J 4 p ) 
(6 1 ,B 1 ),...,(6,,B,) 



is the FOX function. Thus, we can find the solution by substituting Eq.(12) 



into Eq.(4), which yields (see fig.l) 



p(x,t) 



I 2 + 8 



- 2+e 



+00 



dx'x'^ 1 p(x')H^'2 



J\2+8 



(2 



I 1 2+6' 

(0,1), 



2+e ' L > 



(13) 

In Fig.l, we show the behavior of the above equation by considering typical values of 7, Af and 9 with p(x) = 
x 1_A ^(5(x). Note that the solution obtained above leads to an anomalous spreading of the initial condition 
due to the presence of the time fractional derivative and the spatial time-dependent diffusion coefficient. 
The present of Fox function in Eq.(5) may be associated to the changes produced on the probability density 
function for a waiting time that, for this case, has a long tailed behavior, in contrast to the usual one. At 
this point, it is interesting to analyze the asymptotic behavior of Eq.(13) [25]. For simplicity, we consider 
p(x) = x l ~ J ^6(x), so p(x,t) = Q{x,t) and the asymptotic behavior of p{x,t) is 

2 + 6 ,_i *i 1 . 1 



p(x,t) 



2T(M/(2 + 9)) 



(2 — 7)~2^(2+8)(2- T ) 2( 



, AT(-y-i) 



(2 + 6) 2 Dr» 



) (2 + 0)(2- T ) j^j 2-t 



x exp{— (2 — 7)72-7 (- 



12+0 



(14) 



{2 + 9) 2 Dt"f 

In this direction, Eq.(14) can be considered as an extension of the asymptotic behavior of homogeneous and 

isotropic random walk models [26]. 

Case 2. D(t) = , i-e. D(s) = Ds~ a . 
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Using the same method as in case 1, we obtain 



G(x,t) 



and 



2r(^) '(2 + 0) 2 Dfr+> 



-)wHf 



2.0 



12+0 



(2 



)2 Dt y+a 



(1 



2+8 ■ 



7 + a) 



(0,1), (^#,1) 



A' 

2-\- 6 



4-oo 



dx'x'^pix') 



(15) 



x H% 



\x-x'\ 



i\2+e 



(2 + 6) 2 Dt->+ a 



(0,1), C 3 ^,!) 



(16) 



Let us go back to Eq.(l), and consider the external force F (x) oc x\x \ a -\D(x,t) = D6{t)\x\- e and 
fx = 2, v = 1. In this case , analytical solution can not easily be obtained for a generic a, 8. However, for 
6 0, and a + 6 + 1 = 0. By following the same procedure as the one in the above case, an exact solution 
can be obtained and it is given by 



p(x,t) 



1 



2T(^±$£)r(2±^) (2 + 9) 2 Df< 



2+8 



\2+8 



(2 + e) 2 Dv 



(l - J±L y) 

< L 2+6" I) 



<(2+e)D' '-h \ 2+8 ' L > 



(17) 



where, for simplicity, we are considering the initial condition p(x, 0) = 5(x), and the external force (drift) 
F(x) = /C.t|x| q_1 . The second moment is given by < x 2 >cc i^+e, which scales with the exponent 
and clearly depends only on 7 and 9. So, when < 1, = 1 and > 1, the system is sub-diffusion, normal 
diffusion and supper-diffusion respectively 

Let us now discuss Eq.(l) in the presence of the external force by employing a mixing between the 
spatial and time fractional derivatives. For simplicity, we consider (Af, v 1 0) = (1, 1, 0) and F(x) = —ICx. For 
this case, Eq.(l) reads 

§j^p(x,t) = j\t'D{t-t')^p{xX) - ^[F(x)p(x,t)}. (18) 

Applying the Fourier and Laplace transform to Eq.(l) and employing the Riez representation for the spatial 
fractional derivatives, we have 



s~<p(k, s) - 8 r-1 £(fc, 0) = -D(s)\W l p(k, s) - JCk—p(k, s), 

dk 



(19) 



where p(k,t) — F{p(x,t)} = p(x,t)e lkx dx , so p(k,0) = 1. Here, we consider the diffusion coefficient 
given by D(t) = D ^ (a ) , i.e. D(s) = Ds~ a and p{x, 0) = S(x). So we have 



p(M = £-( 

C / J7.I 



1 r D(s)\k\\ n W s-r- 1 

~~!\ ) e 1^- 

n\ JCp s 7 + nfifC 



(20) 



In order to perform the inverse of Laplace transform on Eq.(20), we express e~ z in terms of Fox function, 
i.e. e~ z = H \ \z I (0,1) ] ■ Then we can get the Laplace inverse of p{x, s) by applying the property of the 
Laplace inverse of Fox function [25] 

rt 1 r D\k\H' a 



1 r 1 



ICfJ, 



(0,-a) 
(n,l) 



(21) 
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where, E a _p(x) is the Mittage-Leffler function defined as E a _p(x) = X^riLo r(na+p) \^\- Here, we used the 
property of the Laplace transform of convolution formula, i.e. C[f*g] = f -g, where f*g = J* dt' f(t — t')g(t'). 
Note that the solution of Eq.(18) is a stationary one given in terms of Levy distributions. This feature is a 
characteristic of the presence of the spatial derivatives in the diffusion equation which changes the probability 
for a jump length (see [27] and references therein). In order to get the solution for Eq.(18), we need to get 
the Fourier inverse of p(k,t). Therefore, we only need to perform the inverse of Fourier transform on 
H 



1.0 

1.1 



' D\k\H t<x 


(0,-cO' 




(n,l) 



By employing the procedure presented in [28] , we have 



and 



D\k\n h 



TV 



\ x \ ( Sl/U 

2 [ Dt' a> 



(- -) a -) 

\2> 2/' V ' 1> 



p(x,t) = 



dt' 



1.1 



n\py^\x\t' H2 - 2 



\ x \ , ^ si/u 
2 { Dt' a> 



V 2 ' 2 /> V ' 2 / 



(22) 



x E a , 1 (-nfj l ]C(t-t'y). 



(23) 



Now, we consider a particular case of Eq.(l) for 7 = 1 and nonzero values of p and 0, and consider a 
linear drift, i.e. F(x) = — Kx. For simplicity, we employ D(t) = DS(t) and the initial condition p(x, 0) = 6(x), 
then Eq.(l) yields to 



m p{x,t) 



1 d 



-W,t)] U ) ^ T ^{x^- 1 F(x)p(x,t)} 



dx^ 



(24) 



Let us investigated time dependent solutions for Eq.(24). We use similarity methods to reduce Eq.(24) 
to ordinary differential equations. The explicit form for these ordinary differential equations depends on 
the boundary conditions or restrictions in the form of conservation laws. In this direction, we restrict our 
analysis to find solution that can be expressed as a scaled function of the type [29] 



p(x,t) = 

Inserting Eq.(25) into Eq.(24), we obtain 

K 4>{t) 2 <t>(t)'dz lZ p[ )l 

By choosing the ansatz 



1 



P(*) 



Fl 



D 



0( t yV(„-i)+fl +|i+1 dz 



JC kD 



(25) 



(26) 



(27) 



0(i) 2 <f>(t) 0(t)^("-l)+9+H-l ' 

where k is an arbitrary constant which can be determined by the normalization condition. By solving 
Eq.(27), we have that 



*(t) = [(#0))«e-*™ + ^(l 



(28) 
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where £ = N{v — 1) + 6 + [i. By substituting Eq.(27) into Eq.(26), we obtain 



dz dz^~ x dz 

Then, we perform an integration and the result is 



'~\p(zr} = -k^p(z)]. (29) 



z 



-1-6 ^ 1 - 

dz^ 



-p{z) v = -kz^p{z)+C, (30) 



where C is another arbitrary constant. Also, we use the following generic result [30]: 

Di[x%a + bxf]=a s ^ [a + l] x a - s {a + bxf- 5 (31) 
1 [a + 1 — o\ 

with D & x = d s /dx s and 8 = a + /3 + l. By defining g{x) = (a + bx)» and A = a(l— ^) — 8, and rearranging 
the indices, Eq.(31) can be rewritten as follows: 

d 5 Mx)Y = aa r[a + i-V Ag(a:) ' (32) 

For this case, we consider the ansatz p(z) = Az^ (1 + bz)». By using the property of Eq.(32) in Eq.(30) and 
,for simplicity, choosing C — 0, we find 

(2-p)(p + 9) 
a \-2p-9 ' 

2-p 

v = 



1 + fX + i 



In this case, we have 

PM = m* [ (l + bz)i^ + » + e^ -°> ^ 

where 0(f) is given above, .4 = [~fc r(a+i) ] 1 ~ 2f '~ 8 an d b is an arbitrary constant (to be taken, later on, as 
±1 according to the specific solutions that are studied). Several regions can be analyzed. For simplicity, we 
illustrate two of them: — oo < p < — 1 — 9 with 9 > 0, and < p < 1/2 with < 9 < 1/2 — p. Let us start by 
considering the region — oo < p < —1 — 9. Without loss of generality, we choose b = — 1. The normalization 
condition implies (see Fig. 2) 

•A b r-m — ^ ] 1 - 2lJ -"dz = 1. (35) 



So 



rri - u - 9} 

-pi- p 2 +^e-26'-2^ 1 -p|- i-^+p 2 +e 2 +2^g - l • y ■> 
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FIG. 2: Behavior of <j)(t) N ' p(x, t) versus x/<f)(t), which illustrates Eq.(32) with typical values for /i and 9 satisfying 
< n < — 1 — 9 and 9 > 0. we notice that the distribution vanishes at the abcissa equal ±1, and remains zero outside 
of this interval. 
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FIG. 3: Behavior of cj>{t) N p{x,t) versus x/<j)(t), which illustrates Eq.(32) with typical values for /i and 9 satisfying 
< n < 1/2 and < 6 < 1/2 — fj,. we notice that the distribution vanishes at the infinity. 



Let us now analyze the region < \i < 1/2 with 0<#<l/2 — /i. Again without the loss of generality, 
we choose b = 1. The normalization condition implies (see Fig. 3) 



^4 



T^ r l+e-At 2 -At6> i 
1 I 1-2^-0 i 



(37) 



Let us finally mention a connection between the results obtained here and the solutions that arise from 
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the optimization of the nonextensive entropy [31]. These distributions do not coincide for arbitrary value of 
x. However, the comparison of the \x\ — > oo asymptotic behaviors enables us to identify the type of Tails. 
By identifying the behavior exhibited in Eq.(34) with the asymptotic behaviors l/\x\ 2 ^ q ~ 1 ^ 1 that appears in 
[31] for the entropic problem, we obtain 

9 = 3 + " + l (38) 

This relation recovers the situation for 6 = 0. 



III. SUMMARY AND CONCLUSIONS 



We have analyzed a generalized fractional diffusion equation which presents the spatial and time frac- 
tional derivatives, includes a linear external force F(x) = —JCx and a spatial time-dependent diffusion 
coefficient D(x,t) — D(t)\x\~ e , also takes A/"-dimensional into account. By using Laplace and Fourier trans- 
form, the Green function method and normalized scaled function, we can find the explicit solutions p(x, t) 
which subjects to the natural boundary condition p(±oo, t) = and the initial condition p(x, t) = p(x). In a 
word, we have extended the results previously obtained by the other authors by inviting an linear external 
force, a spatial time-dependent diffusion coefficient and A/"-dimensional case. We have also discussed the 
connection with nonextensive statistics, providing the relation between our solutions and those obtained 
within the maximum entropy principle by using the Tsallis entropy. Finally, we expect that the results 
obtained here may be useful to the discission of the anomalous diffusion systems where fractional diffusion 
equations play an important role. 
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